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1. INTRODUCTION 


The basic objective of the research project was to determine the mass 
required to protect a surface during atmospheric re-entry by injecting a 
cold gas tangentially along the surface. Such a thermal protection system 
is known as film or slot cooling and is of considerable utility in many 
areas. It involves the mixing of the high enthalpy external stream with 
the coolant in the vicinity of a solid wall which is the surface being 
protected. A boundary layer forms along the wall, and eventually inter- 
acts strongly with the external mixing zone. 

A consideration of this problem revealed that the boundary layer 
equations could be employed to approximate the most important flow con- 
figurations which occur in film cooling. However, even with this simpli- 
fication, the equations to be solved are complex, non-linear, coupled 
partial differential equations. Due to the relative weakness of classical 
methods of boundary layer analysis, finite-difference methods were chosen 
as the best approach to the formidable problem under consideration. 

A literature search revealed that the use of finite-difference tech- 
niques for the solution of the boundary layer equations has become popular 
only recently. Discrepancies and inconsistencies are common among recent 
publications. Many of the proposed techniques required an initial condition 
in the normal component of the velocity, v. This requirement is inconsistent 
with the boundary layer approximations. Other methods employed transformations 
which, in many cases, restricted the range of applicability of the technique. 

Instead of blindly adopting a proposed technique and trying to force it 
to our problem, a more systematic study was deemed to be not only desirable 
but necessary. Otherwise numerous dead ends would probably have been en- 
countered, and if a dumber 1 were obtained, it could have contained cata- 
strophic errors. Also, the 'number 1 would have been the only product of the 
extensive efforts and the large digital computational costs. Necessary' 
simplifications and changes in interest which often occur, could have meant 
that even correct answers would have been of little value. Hence, a rela- 
tively basic study of the application of finite-difference methods was under- 
taken with the objective of developing an algorithm capable of obtaining 
accurate solutions to film cooling problems. 



An earlier report. Reference 1, contains the details of the proposed 
finite-difference algorithms and how they differ from those given in the 
literature. Numerous solutions were presented and graphical and tabular 
data were provided to establish the accuracy of the procedures. Further 
extensions including applications to turbulent flows are contained in 
Section 3 of this report. 

As a consequence of the systematic study which was undertaken, several 
significant advancements in the general theory of finite -difference methods 
were made. The main items are: 

1. Procedures for generalizing finite -difference solutions of 
parabolic partial differential equations, 

2. Criteria for choosing a priori suitable increments in the in- 
dependent variables for effecting accurate solutions to such 
equations , and 

3. Procedures for reducing the computational effort beyond that 
accomplished through Item 1 (reductions of several orders of 
magnitude are common). 

This material is presented in Section 2. In general, the procedures and 
criteria are not limited to any specific algorithm or problem; therefore, 
they should be of value to most investigators engaged in finite-difference 
solutions of partial differential equations of the parabolic type. 

Section 4 contains results from the study of film cooling problems. 
Numerous results for incompressible and compressible laminar film cooling 
have been obtained. The initial turbulent results which have been ob- 
tained are not presented in Section 4. It was felt that it might be mis- 
leading to present the meager preliminary results currently available with- 
out further refinement. The investigation has not been completed; several 
assumptions need to be examined more carefully. The study will be con- 
tinued, however, because it represents the subject of Mr. P. L. Kong's 
Ph.D, dissertation which is still incomplete. Additional publications con- 
taining further details and extensions are also planned. 
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2. GENERAL ADVANCEMENTS IN FINITE DIFFERENCE METHODS APPLICABLE 
TO PARABOLIC PARTIAL DIFFERENTIAL EQUATIONS 

2.1 Similarity and the Generalization of Finite-Difference Solutions 

Reference [2]* 5 a publication by one of the authors 9 contains a de- 
tailed discussion of the subject of similarity and the generalization of 
finite difference solutions. This section is a synopsis of [2], and it 
provides the necessary equations for the analyses and discussions given 
in the latter portions of the report . 

Certain sets of partial differential equations can be reduced to 
sets of ordinary differential equations by appropriate transformations 
of the independent and dependent variables. Such solutions are generally 
designated as "similarity solutions" and are called "similar" because 
the profiles can be made congruent if they are plotted in coordinates 
that are made dimensionless with reference to the appropriate scale fac- 
tors . 

The search for similarity transformations and techniques of solv- 
ing the resulting ordinary differential equations has been extensive. 
Impetus to the search for similarity transformations has been provided 
mainly by the supposed mathematical simplifications. For example 9 
Schlichting [4] states that the reduction of a system of partial differ- 
ential equations to one involving ordinary differential equations "evi- 
dently constitutes a considerable mathematical simplification of the 
problem." However 9 the ordinary differential equations that result are 
usually non-linear and effecting their solutions is difficult. 

It was shown in [2] that the simplifications and advantages of simi- 
lar problems are generally not lost if the solution to the original set 
of partial differential equations is effected with finite-difference 
methods in the physical plane . Furthermore 9 the similarity transforma- 
tions can often be deduced by examining the governing difference equations 
written in terms of a generalized difference operator. If this property 
exists and is recognized 9 a large savings in computational effort is 
usually realized. 

Since few practical problems are similar 9 the most significant 


‘“Numbers in brackets refer to entries in Section 5 9 REFERENCES. 
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findings was that considerable savings can also be effected in non- 
similar problems . Not just one solution , but a whole family of solu- 
tions can often be effected with a single integration in the physical 
plane* One of the key steps in the generalization is to set up the 
problem such that the numerical calculation can be made without specify- 
ing the increments in the independent variables* Although the solution 
is effected in the physical plane, n ew dependent variables must some- 
times be introduced in order to accomplish this objective* For example, 
the quantity v 1 , defined as vAx/Ay, might be used in the numerical in- 
tegration in place of the dependent variable v. Such quantities lead 
to the generalized dependent variables (the so-called similarity variables 
in similar problems) later in the analysis when the unspecified incre- 
ments 9 Ax and Ay, are eliminated* 

A parameter M that specifies the ratio between the spatial incre- 
ment squared h 2 and the time increment At arises in all finite -difference - 
algorithms of parabolic equations because of the nature of such equations* 
The definition of the particular parameter used in the analysis is de- 
pendent on the problem because other constants are introduced in order to 
reduce the parameter to a dimensionless quantity. The parameter is, of 
course, unique to the numerical solution. A key in the present analysis 
is the fact that the solution obtained with a stable and consistent 
finite-difference approximation in the limit as h and At approach zero is 
independent of h 2 /At; that is, h 2 /At influences the results obtained with 
such finite -difference approximations, with finite values of h and At, 
only because this quantity influences the error. 

The procedure for effecting the generalization of finite-difference 
solutions can be summarized as follows: 

1. Set up the problem without specifying the increments in the 
independent variables . 

2. Calculate N steps where N can be chosen arbitrarily. 

3. Eliminate the increments from the dependent variables by 
introducing M, N, t (or the analogous spatial variable) and 
possibly other parameters. 

4. For heuristic purposes, separate the solution into two parts; 
the exact solution which is independent of M, and the truncation 
error which is dependent on M. 
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5 . 


Relate the indexes i, N and possibly others to the independent 
variables . 

60 Reduce the number of variables by utilizing the fact that 
only the truncation error is dependent on M. 

The mechanics become much clearer if one examines an actual prob- 
lem 0 One of the simplest parabolic equations will be considered next; 
results from many additional examples are given in [2] and [3]. An 
infinite plate, which is initially isothermal and which suddenly be- 
gins to exchange heat by convection at one boundary, is the problem 
of interest. This is a nonsimilar problem, and one of the limiting 
cases of this problem, which will also be considered, is similar. The 
governing equation is 


1_ 3T _ 3 2 T 
a 3t » 2 

9 y 


( 1 ) 


where the constant a is the thermal diffusivity and T is the tempera- 
ture. The initial and boundary conditions are 

T = T at t = 0 5 0 < y < L 

o ’ J — 


-k 


3T 

3y 


y-Q 


= h(T - T), t > 0 

s 


9T 

9y 


= o. 


t > 0 


y=L 


> 


j 


( 2 ) 


In replacing the derivatives in Eqs. (l) and (2)*by difference quotients, 
the difference equations will be written in terms of a difference oper- 
ator V. V represents any conceivable operator and is introduced in order 
to consider simultaneously all possible finite -difference representations 
of Eq. (3). The operator is not unique. If the discrete variables are 
defiped as 


t = n(At) 

n 

y. = (i - l)Ay 
T. = T at y. and t 

i , n J i n 

the finite-difference representation of Eq. (1) becomes 
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v x* V 2 T* 

1 t i , n y i , n 


where 


V T* = 4 V 2 T* 

t i , n My i , n 


M = and T* = =- 

otAt T 


The initial and boundary conditions in discrete form are 


V x* 

y 1, n 


v T* 
y P. n 


(T* - 1) = Bi+CT* n 

1 . n X , n 


where L = (p - l)Ay and Bi 1 " = hAy/k. Equations (3) and (4) reveal that 
T* n = f 1 (i 3 n 3 M 3 p 3 Bit ) (5) 

Equation (5) shows that the solution can be effected without specify- 
ing either Ay or At although quantities such as p 5 and Bi must be 
chosen before the numerical solution can be effected* Leaving Ay and 
At unspecified is 9 of course , unnecessary but it is a distinct advantage 
in some problems and is the key to eliminating redundant calculations* 
Equation (5) represents the results from Step 1 of the numbered pro- 
cedure* 


If one proceeds N-time steps (see Step 2), one obtains 
T* N = fjU.N.M.p.BiO 


Since this problem is relatively simple 9 the dependent variable does 
not contain any unspecified increments; hence, Step 3 is unnecessary. 
Otherwise it would lead to the relevant dependent variables* 

Next, the solution is separated into two parts: the exact solution 
which is independent of M, and the truncation error, E , which is de- 
pendent on M (see Step 4). 
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T* N = f a (i,N,p,Bit) + E T (i,N,M,p,Bi t ) 


( 6 ) 


The numerical integration yields f ; f is introduced for heuristic 
purposes. It is assumed throughout Section 2 that gross and round-off 
errors are negligible. 

Since Ay and At are unspecified, the meaning of T* in the y-t 
plane remains to be determined (see Step 5). The relationship between 
the index i and the independent variables y and t follows from the 
definition of y = (i - l)Ay if M and t N /N are introduced to eliminate 
Ay and At, respectively. The result is 



or 


n : 


i, N 


Likewise , 





n 


P, N 



(p - 1) 



(7) 


( 8 ) 


If Ay is eliminated from Bi* by introducing p, the Biot modulus follows 


Bi = ~ = Bi t (p - 1) 


(9) 


Alternatively, Ay could have been eliminated by introducing *\[ M/N 1 y 


at. 


Bi' 


h VatT _ g. t J N 



( 10 ) 


If n , ri and Bi are introduced as new variables in place of N, d 

i t n 5 p, N r 9 r 


and Bi* , respectively, Eq. (6) becomes: 


T |% * f 3 <i > n i,N> rl P , N- B1) * E T 


(ll) 


Since f is independent of M, M can be chosen for heuristic purposes to 
be proportional to N. Equations (7) and (8) shows that r\. and T) 

r 1 r N p p N. 

are independent of N for this case; hence, Eq. (11) becomes (see Step 6) 

t* = f 4 (n. ,n p ,Bi) + e t (12) 
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Alternatively, one could have used Bi' and obtained 

T* = f 4 (n. 5 n p 3 Bi f ) + E t (13) 

Because N and t N are free to take on any value, the subscript N is de- 
leted from T* t , ri. and J] It is to be noted that the new 

1 , N 3 N 1 f N p, N 

independent variable, rj ( =y / V at 1 ) , arose from the governing equation, 
and it would have arisen regardless of the geometry, the initial con- x 
dition, or the boundary conditions. This variable will be called the 
diffusion variable and is a meaningful dimensionless coordinate which 
governs the diffusion process. It is identical to the similarity variable 
in similar problems; however, its use and significance is by no means 
limited to similar problems. Equation (12) does not represent a unique 
result. A general discussion of alternative procedures is given in [2]. 

The boundary conditions and the finite length of the region L re- 
sulted in the parameters r\ and Bi. If L is large relative to / ycct , 
physical considerations show that the solution becomes independent of 

D , and one obtains from Eq. (13) 
p 

T* = f(n. ,Bi f ) + E t (14) 

Thus it is immediately clear that the whole family of solutions for this 
limiting case can be obtained with a single numerical integration in the 
physical plane. If h is large relative to k/V at, the solution becomes 
independent of Bi 1 ; that is, for a semi-infinite solid with a step 
change in surface temperature , the solution is 

T* = f(n. ) + E (15) 

l i T 

Tills limiting solution is obviously a similar solution. If the finite 
size of the region is of importance but Bi is large, Eq. (8) becomes 

T * = f (nj ,n p ) + e t (16) 

In the boundary layer examples given in [2], the diffusion variable, 
q, which resulted had an identical definition in the discrete plane 



q and M are of course different, they are: 
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n = y 



and M 


uJAy) 2 

vAx 


The definition of was dependent on the particular problem. 

The conclusions given in [2] are: 

1. Similarity variables , if such quantities exist, can fre- 
quently be deduced from an analysis of the discrete form 
of the original differential equations ahd the respective 
boundary and initial conditions. 

2. The generalization of the governing finite-difference 
equations can save considerable computational effort 
through the elimination of redundant calculations. The 
proposed analysis is applicable to both similar and non- 
similar problems. 

3. Many similar solutions can be effected easily by solving 
the original partial differential equations in the physical 
plane by finite-difference methods. If numerical pro- 
cedures are to be employed in effecting the solution, the 
supposed mathematical simplification introduced with the 
application of similarity transformations must be seriously 
questioned. 


2.2 Practical Criteria for Choosing, a priori , the Increments in the 
Independent Variables 

A major deficiency of numerical techniques is the inability to 
choose increments in the independent variables which will yield a re- 
liable result. Young [5], in his survey of numerical analysis of el- 
liptic and parabolic partial-differential equations, simply stated: 

"None of the papers which have been written on error estimation [35, 
47-49,51,52,65] can be said to come close to satisfying the needs of 
the practical computer user." A new criterion was proposed in [3] 
which partially alleviated this problem and many numerical results 
were presented in order* to establish the validity of the proposed 
criterion. 

The criterion essentially provides a practical means of taking 
into account the magnitudes of the derivations in choosing the increments 
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in the independent variables. A summary of the basis -for the criterion 
follows and provides a starting point for the extensions given in 
Section 2.3, The reader is referred to [3] for additional information. 

Consider the diffusion of heat, mass, or momentum into an in- 
finite or semi-infinite medium. If the diffusion layer is thick, the 
spatial derivatives obviously are small. In addition, the mechanism 

\ 

of the diffusion process causes the time* derivatives to be small when- 
ever the spatial derivatives are small. Hence, large spatial and time 
increments will give completely satisfactory results if the solution is 
desired at relatively large times. On the other hand, if the diffusion 
layer is thin, the spatial and time derivatives are large and the spatial 
increments in the region of influence and the time increments might have 
to be several orders of magnitude smaller in order to effect a meaningful 
solution. The key to a practical error criterion is that the approximate 
magnitude of the derivatives can be deduced from the thickness of the 
diffusion layer. The fact that the thickness of the diffusion layer is 
unknown is surmounted by relating it to a thickness in terms of a new 
dimensionless independent variable, the diffusion variable. Although the 
dimensional thickness of the diffusion layer may vary drastically with 
respect to time or the independent variable analogous to time, the thick- 
ness of the diffusion layer in terms of the diffusion variable, r| , is a 
constant for a similar problem and usually varies over a relatively small 
range in non-similar problems* In addition, the diffusion thickness does 
not vary drastically from problem to problem. The diffusion thickness 
(see [3] for a precise definition) only varied between 1.4 and 5 for the 
diverse examples considered in Table 1 of [3] which included heat con- 
duction problems, natural-convection boundary layer flows, and incompres- 
sible and compressible hydrodynamic boundary layers on flat plates and 
wedges. A notable exception is thermal boundary layers in liquid metals. 

The hypothesis on which the error criterion is based is: The ac- 

curacy of a numerical solution can be estimated from the magnitude of the 
increment in the diffusion variable Ar| , and this parameter provides a 


^Throughout the discussion, one of the independent variables is assumed to 
be time. This is not always the case; e.g. , consider two-dimensional 
steady-state boundary- layer problems. 
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means of intelligently choosing spatial and time increments for efficiently 
effecting finite -difference solutions of parabolic partial-differential 
equations in the physical plane. The number of grid points in the dif- 
fusion layer is the diffusion thickness divided by Arj » hence, a limit on 
the maximum value of Ari fixes the minimum number of grid points which lie 
in the diffusion layer. It should be stressed that the solution is ef- 
fected in the physical plane. The diffusion variable and the correspond- 
ing diffusion thickness are introduced in order to establish an error 
criterion. 

The implications of a restriction on the magnitude of Ari only be- 
come meaningful after a relationship is established between this quantity 
and the parameters of the discrete, physical plane. The procedure for 
determining this relationship was given in Section 2.1. Although the 
cases presented in [3] covered a wide range of physical problems, the re- 
lationship between the numerical parameters and the diffusion variable X] 
was identical in all cases and is 



Hence, a criterion which restricts the magnitude of Arj,and in this man- 
ner guarantees several grid points within the diffusion layer, restricts 
the magnitude of V^/n' . The proposed criterion is 

(16) 

Criterion (16) guarantees between 5 and 18 grid points in the diffusion 
layer at the location of interest, N(At) or N(Ax), for the set of prob- 
lems given in Table 1 of [3]. The criterion is designed to limit the 
truncation error and is not applicable if other errors are of importance. 
Round-off errors, programming errors, etc., can be accessed easily in 
many Engineering problems by applying overall balances to the system, 
e,g., an energy balance, a force balance, or a mass balance. 

Criterion (16) provides a practical means of choosing the incre- 
ments in the independent variables for the numerical integration. It 
is simple to use and is applicable to a vast number of problems in a 
variety of fields. In the majority of cases investigated, the numerical 
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percent error lies between 0.2 and 2 percent at Aq = 0,3. The main ex- 
ceptions, were quantities which were asymptotically approaching zero; 
obviously 3 the percent error becomes meaningless in such cases. 

If M/N were the only influences on the accuracy of the finite- 
difference solutions 3 the values of the increments should be chosen 
such that small values of M are obtained. For example, M = 4 and 

\ 

N = 25 would result in the same value of Ari as M = 40 and N = 250 , but 
the latter case would require an order of magnitude more computational 
e f f ort . 


Expressions for truncation errors would indicate that one should 
limit the increments in the independent variables Ay and At. In 
contrast , the present error hypothesis places limits on the sizes of Aq 
and M. Aq was introduced in order to take into account the magnitude of 
the derivatives; the minimum value of M must be restricted in order to 
avoid violations of physical laws. Criteria have been derived to avoid 
potential violations of physical laws for a variety of finite-difference 
methods [6,1]. These criteria are of the form 


M > 

for heat-conduction problems, and 


, N > S 2 


(17) 


(18) 


for boundary -layer problems. An additional criterion arises in boundary- 
layer problems if the energy equation is simultaneously being considered. 
The value of £ is a function of the finite-difference method and the prob 
lem. It is of the order of one for problems with step changes in the 
neighborhood of the change. If physical laws are being violated, meaning 
less oscillations will occur. These can easily be detected from the re- 
sults as long as they are monitored such that successive samples include 
both odd and even values of n. Hence, the form of the criterion and the 
order of magnitude of £ is suffienct for the practical computer user. 
Criterion (18) shows that the numerical value of £2 is of limited value 
anyway because u* ^ is normally unknown; this is a characteristic of non- 
linear problems. 


Numerous results were provided in [3] in order to establish the 
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utility of Criterion (16). These results will not be repeated here. 

The conclusions drawn in [3] from the numerous and varied test cases 
are : 

1. The truncation error is strongly dependent on An and rela- 
tively independent of most other parameters. 

2. If M is sufficiently large to avoid violations of physical 
laws, it has a minor influence on the accuracy of finite- 
difference solutions effected with the Crank -Ni cols on 
algorithm. 

3. The criterion provides a priori means of choosing the in- 
crements in the independent variable which are required in 
order to effect a meaningful solution. The increments are 
directly related to the time or x- location of interest and 
to the relevant molecular diffusivity. 

4. When the criterion was satisfied , engineering accuracy was 
obtained in all cases considered with the possible exception 
of compressible boundary layers at high Mach numbers. The 
reason for this exception was given. 

5. The criterion gives satisfactory results for non-similar 
problems , for non-linear problems, and for diffusion in 
finite regions. It should be used with caution, however, 
in areas which are vastly different from the examples used 
in the present verifications. Although the form of the 
criterion should remain applicable, the constant 0.3 may prove 
to be too small or too large. 

2.3 The Use of Nonuniform Grids in Finite-Difference Sol uti ons— Numeri cal 
Integration with Constant Accuracy 

The accuracy of a finite-difference solution to a diffusion problem 
usually increases as time* increases if the increments are held constant. 
Usually the derivatives are becoming smaller; hence, expressions for the 
local truncation error clearly show why the error decreases in such cases. 


*Again , throughout the discussion in Section 2.3, one of the independent 
variables is assumed to be time. This is not always the case; e.g., x is 
the analogous independent variable in the boundary layer problems. 
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Small spatial and time increments are generally required in order to 
effect an accurate solution at small times, whereas large spatial and 
time increments give completely satisfactory results if large times are 
of interest. 

The procedure given in Section 2.1, which enables one to leave 
the increments unspecified, greatly allevaiates this difficulty in 
some cases by leaving the increments unspecified. Obviously, the re-' 
suits obtained after N steps can be applied at any time if the time incre- 
ment has not been specified. Although the magnitude of the problem has 
been reduced by leaving the increments unspecified, not all difficulties 
and expensive calculations have been eliminated. One area which leads 
to costly calculations is nonsimilar problems if large ranges of the 
additional parameters are of interest. For example, consider the 
parameters ri and Bi in Eqs. (12 and (14), respectively, of Section 2.1. 
Even though a whole family of solutions is obtained with a single in- 
tegration, the calculation still can be expensive. A turbulent boundary 
layer flow presents difficulties because small increments are required 
in order to include grid points in the non -turbulent sublayer. However, 
the thickness of this layer is small relative to the thickness of the 
boundary layer. The use of a uniform network is clearly unacceptable 
in this case. 


Furthermore, leaving the increments unspecified in no way removes 
the dependency of the accuracy on the number of time steps, N. Ideally, 
a procedure which results in a solution whose accuracy is independent 
of the number of time steps is desired. 


The use of a non uniform grid is by no means a new idea. Many 
investigators simply double or triple the time and/or spatial increments 
every so often. For example, Pletcher [7] used the following distri- 
bution in the y-increment in his finite-difference solution of a turbu- 
lent boundary layer. 


Ay + = 4 


0 < y 


80 


Ay 


10 


Ay = 20 
Ay *" = 40 


Ay = 100 1180 < y 


80 < y _< 180 

180 < y + _< 380 

380 < y + _< 1180 

+ 
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A procedure like Pletcher’s suffers from two main disadvantages . 

First 3 the changes in the local truncation error are large and abrupt. 
Since the local truncation error due to the y-derivatives is of the 
order of (Ay) 2 in Pletcher’s algorithm., a change by a factor of 
either 4 or 6.25 occurs from one region to the next. Second 3 the local 
truncation error at the interface between any two regions is, in general , 
larger than that in either region. Hence 9 some of the desired gain is 
offset by the large truncation error which occurs at these interfaces. 
Abrupt changes in the time increments usually suffer from only the first 
of these two disadvantages. This is because most numerical procedures 
employ forward difference quotients for the time derivatives. 

Smith and Cebeci [8] appear to have been one of the first to use 
a more smoothly varying grid. They employed a geometrical progression 
for the variation of the increment in their transformed coordinate nor- 
mal to the surface in a turbulent boundary layer study. Calling this 
coordinate y 5 successive y increments were 

Ay x , rAy x , r 2 Ay x , r 3 Ay, , etc. 

Hen ce 5 


Ay. = rAy. _ x , i > 1 
or 

Ay. = r I-1 Ay 1 , i > 1 
where r _> 1. It follows that 

y = Ay 1 + rAy x t r 2 Ay L + ••• + r 1 1 Ay 1 


or 


(19) 


( 20 ) 



Smith and Cebeci did not discuss the influence of r on their accuracy 5 
nor did they indicate any basis for choosing r or a geometrical pro- 
gression other than the following: "This variable grid system permits 

shorter steps close to the wall and larger steps away from the wall, 
and thus maintains computing accuracy." They employed an arbitrary 
spacing in the other transformed coordinate. More recent studies by 
Harris [9] and Adams [10] also used a geometrical progression and 
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credited Smith and Cebeci for proposing it. 


It is interesting to note that so far no one has been successful 
in providing a rational, systematic means of varying the grid spacing. 

A basis or criterion for establishing the manner in which the incre- 
ments should be varied has not been found. 

The objective is to vary the spatial and time increments during N 
the integration such that the accuracy of the results is independent 
of both the time and the spatial location. Implicit in this objective 
is the requirement that the process remain stable and satisfy the funda- 
mental laws governing the phenomenon 9 e.g. , the second law of thermo- 
dynamics. It is probably impossible to fully accomplish this objective. 
However, the criterion recently developed by Clausing [3] appears to 
provide a means of approaching this subject in a rational manner. Hereto- 
fore, a basic ingredient, a practical error criterion, was missing. 

It was established in [3] (also see Section 2,2) that the accuracy of 
finite- difference solutions is strongly dependent on the increment in 
the diffusion variable. An. The proposed criterion is: 



The relationship between the diffusion variable and the parameters of 
the discrete physical plane can be established with the procedures out- 
lined in Sections 2.1 and 2.2. The stated objective when translated into 
Criterion (16) through Criterion (18) requires that the calculation pro- 
ceed in a manner such that Aq remains relatively constant. Simultaneously, 
one cannot allow M to become too small. Clearly, the y-increments must 
be dependent on both y and time. Also, the time increment cannot be 
chosen independently, but is strongly dependent on the variation of Ay } 
with time. Smith and Cebeci [8], Harris [9], and Adams [10] use incre- 
ments in the transformed variable normal to the wall which were independent 
of the other coordinate. Furthermore, all three of these investigators 
used an "arbitrary" or constant spacing of the grid points in the other 
direction. The two increments were in no way related to each other. 

Since large changes in Ay and At at arbitrary values of i and/or n 
are to be avoided if possible, geometrical progressions in both increments 



will be considered. Geometrical progressions appeal^ to be a natural 
choice for diffusion problems of the boundary layer type. Initially 
Ay will be assumed to be independent of time ; the problem will be 
generalized further later in this section. Hence, 


i -1 


= r Ayj 


and 


a - n- 1 A 

At = r At, , 

n t 1 3 


therefore , 


y 5 = a Yi 


f r 1 - 1 

_y 

r - 1 

y 


and 


i > 1 


n > 1 


( 22 ) 


(23) 


i > 1 


(24) 


n 

r, - i 

t = At I =- 

n li r. - 1 


n > 0 


(25) 


If the "edge" of the diffusion layer lies at i = I and the calculation 
is to proceed to n = N, the equivalent values required with a uniform 
grid with increments of Ay, and At x are, respectively. 


i -l 


I. = 1 +l 


r - 1 

y 


(26) 


and 


/ n 

(r - 1) 

N = 7 ^ 

e (r - 1) 


(27) 


The savings in computational effort which results from these geometri- 
cal progressions can be deduced from Table 2,1. This table shows the 
influence of r and N on /N or r and (I - 1) on (I e - 1 ) / ( I - 1). Typically* 
one would probably use a value of I which is less than 100 and N which 
is less than 1,000; hence, the upper half of the table is designed for 
determining I and the lower half for determining N g . 

Two problems will be used to show the influences of r y and r t on 
the accuracy: heat flow in a semi-infinite solid with a step change in 
surface temperature and incompressible laminar boundary layer flow over 
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Table 2.1 The Influence of r on the Effective Values of I and N 


r 

II \ 

or ( I - 1 j\ 

' 


N/N or (I - 
e e 

1 )/(I - 1 ) 




1.001 

1.002 

1.005 

1.010 

1.020 

1.050 

1.100 

1.200 

5 

1.0 

1.0 

1.0 

1.0 

1.0 

1 - 1 . 

1.2 

1.5 

10 

1.0 

1.0 

1.0 

1.0 

1.1 

1.3 

1.6 

2.6 

15 

1.0 

1.0 

1.0 

1.1 

1.1 

1.4 

2.1 

4.8 

20 

1.0 

1.0 

1.0 

1=1 

1 . 2 

1.6 

2.9 

9.3 

25 

1.0 

1.0 

1.1 

1.1 

1.3 

1.9 

3.9 

18.9 

30 

1.0 

1,0 

1.1 

1.2 

1.4 

2.2 

5.5 

39.4 

35 

1.0 

1.0 

1.1 

1.2 

1.4 

2.6 

7,7 

84.2 

40 

loO 

1.0 

1.1 

1.2 

1.5 

3.0 

n.i 

183 

45 

1.6 

1.0 

1.1 

1 . 3 

1.6 

3.6 

16.0 

406 

50 

1.0 

1.1 

1.1 

1.3 

1.7 

4.2 

23.3 

910 

55 

1.0 

1.1 

1.1 

1.3 

1 0 8 

5.0 

34,2 

2059 

60 

1.0 

1.1 

1.2 

1.4 

1.9 

5.9 

50.6 

4696 

70 

1.0 

1.1 

1.2 

1.4 

2.1 

8.4 

112 

24920 

80 

1.0 

1.1 

1.2 

1.5 

2.4 

12.1 

256 

> 10 s 

100 

1.0 

1.1 

1 . 3 

1.7 

3.1 

26.1 

1378 


150 

1.1 

1.2 

1.5 

2.3 

6.2 

200 

> 10 s 


200 

1.1 

1.2 

1.7 

3.2 

12.9 

1729 



250 

1.1 

1.3 

. 

2.0 

4.4 

28.1 

15864 



300 

1.2 

1.4 

2.3 

6.3 

63.2 

> 10 5 



350 

1.2 

1.4 

2.7 

9.0 

146 


1 


400 

1.2 

1.5 

3.2 

13,1 

344 




500 

1.3 

1.7 

4.4 

28.8 

1996 


i 


600 

1.4 

1.9 

6.3 

65.1 

12051 


! 


700 

1.4 

2.2 

9.1 

151 

74837 




800 

1.5 

2.5 

13.3 

358 

> 10 s 




900 

1.6 

2.8 

19.6 

861 



i 


1000 

1,7 

3.2 

29.1 

2096 



; 

i 
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a flat plate, the Blasius problem. (See examples 1 and 3 of [2] for 
governing equations and other details.) If one chooses M = 10 and 
were to integrate to An = 0.05, a uniform grid would require I - 140 
and N = 4*000. A total of approximately one -half million grid points 
would result if a constant value of I is employed. This could be re- 
duced by approximately one -third if I were varied as the thickness of 
the diffusion layer increased. Using r y = r^ = r^ = 1.05 requires 
I = 45 and N = 109, a total of approximately 4,900 grid points with a 
rectangular grid. This represents a reduction in the computational 
effort of three orders of magnitude. Integrating an equivalent of 4,000 
steps is not unrealistic. If Arj must be less than 0 . 3 in order to 
obtain the desired accuracy, N g = 4,000 would result in a useful solution 
over a time period from t Q to 36 t Q or a Reynolds number ranging from 
Re^ to 36 Re Q . If I = 150 and N = 4,000 were actually used, round-off 
error would probably strongly influence the results. 

The influence of r^ on the accuracy is shown in Figs. 2.1 and 2.2. 
Figure 2.1 shows the variation of the percent error in the heat flux 
through the surface with Arj for the heat conduction problem. Six differ- 
ent values of r y are used. Figure 2.2 gives results for the Blasius 
problem. Both sets of results were obtained with r t (or r ) equal to 
one. Figure 2.1 shows that, relative to the influence of Arj, the in- 
fluence of r^ is small if r^ is less than 1.3. The apparent influence 
of r y is greater for the Blasius problem but the original error curve is 
misleadingly small because it parallels the abscissa. The approximate 
number of grid points required in the y-direction at An = 0.15 for 
r y = 1, 1.05. and 1.3 is 48, 26 , and 11, respectively. 

Figures 2.3 and 2.4 show the influence of r (or r ) for these same 

t X 

two problems. A value of r = 1.05 was chosen in both cases because the 
integration was carried to rather small values of Arj. The integration to 
An = v 0.05 required 1,600 , 91, and 32 steps for r (or r ) equal to 1.0, 
1.05, and 1.2, respectively. The computational effort required to inte- 
grate to Aq = 0.05 with r^ = 1.05 and r. =1.2 was less than 1 percent of 
that required with a uniform grid. The influence of r or r^ on the ac- 
curacy is smaller than the influence of r , provided that the calculation 


19 




20 


Figure 2.1 The Influence of r on Accuracy — Heat Conduction in Semi -In finite Solid 







avoids violations of physical laws. The influence of • r or r is es- 

t X 

sentially negligible compared to the influence of An if r is less than 
approximately 1.1. 

If a nonuniform grid is employed, expressions for ri. and An must 
again be derived. The following definitions will be used: 

(A yi ) 2 

M = (28) 

ctAtj 


(Ay x ) 2 

M E 

i, N oAt 

N 


and 


(29) 


Al1 l E \ N 


1 1 , N 


(30) 


Since Ay j and At J are unspecified, the meaning of T* ^ in the y-t plane 
remains to be determined. The procedure with a nonuniform grid is identi- 
cal to that used with a uniform grid* One starts with the definition 


of y. 9 


Eq. (24) in this case. 


M is introduced in order to eliminate 


Ay x . This introduces 
Eq, (25), This gives 


At which is eliminated by introducing N and t 


y. 


i -1 


n. 


i » N 




r - l i 
y 


y 2 


N 


- 1 


, r - 1 
u t 


(31) 


or if N is introduced 

e 



The definition of Ap^ , Eq. (30) combined with Eq. (32) gives: 



(32) 


(33) 


The only differences between these results and the expressions ob- 
tained for a compressible or incompressible boundary layer on a flat 
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plate lie in the definitions of T) and M . Although Ar^ is the smallest 
increment in this variable, it is the logical choice for use in the 
error criterion, Criterion (16). This is demonstrated by the error 
curves given in Figs. 2.1 through 2.4. M is the smallest value of 
M and should be used in Criteria (17) and (18). 

Equation (33) shows that Ari is by no means constant. Indeed, 
smaller values of Ari are generated much more quickly if r > 1. In 
addition, Eqs. (23) and (24) give 



(34) 


Hence, N is growing smaller with increasing N. Clearly, if Ayj is 
independent of time, our objective cannot be accomplished. Dangerously 
small values of M will be generated and the calculation will become 
inefficient due to unnecessarily small values of Ay . This will, in 
turn, result in such large values of I that round-off error may become 
of importance. 

An obvious choice for meeting the objective is to assume Ay is 
given by 


Ay, - 


^i.n 


The corresponding expression for Ari becomes 



Since 


( 35 ) 


(36) 
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If this value of r is substituted into Eq, (36), Ar^ becomes 

f - 1) H, 1 

An, 


N *1 
r - 1 

t 


( 38 ) 


Equation (38) shows that Ar^ would be strongly dependent on N only if 
N were small. The asymptotic value of Ar^ follows from Eq. (38) and is 


lim An 
N ■+ 00 



(39) 


Thus 9 the asymptotic value of.Ar^ has been changed from zero to a finite 
constant by introducing r and relating it to r^ such that M is a con- 
stant. M 1 and r can be chosen in order to give the desired asymptotic 
value of An 9 i.e, 9 the desired accuracy. 

Although both objectives would appear to have been accomplished 
(at least to some degree) by this procedure 9 the use of r causes a 
serious problem. A nonorthogonal set of grid lines results. The grid 
is orthogonal only at i = 1 and becomes more and more skewed with in- 
creasing i. If this skewness is ignored 9 serious error can result in 
some problems. To correctly take the skewness into account would be 
extremely difficult or impossible for complex sets of equations. Thus 9 
the use of r^^ appears to be impractical. 

Since any continuous variation of Ay with N would result in a 
skewed network 9 a periodic abrupt change appears necessary in order to 
approach the main objective. The proposed procedure is to periodically 
drop every other grid point in the y-direction such that the value of 
Ari l will not become less than some prescribed value 9 An c . When every 
other grid point is dropped 9 An w ill approximately double; hence 9 the 
accuracy will not be constant but it will only vary over a narrow range. 
Criterion (16) suggests An c should be approximately 0.1 or 0.15. 

If the superscript m is introduced to identify the successive 
values of Ay i 3 and Ay* E Ay^ 9 then: 

A m > „ m-1 x , m-1 , N 

Ay = (1 + r ) Ay, , m > 1 (40) 

1 y 1 
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Deleting every other grid point results in another geomtric progres- 


sion but with a different value of r 


Successive values of r are 

y 


m , m-1 n 2 

r = (r ) 
y y 


m > 1 


(41) 


r^ must be chosen with this progression in mind. In general, the calcu- 
lation should probably be limited to three or four changes. Equation (40) 
shows that M will change by a factor of (1 + r™ 1 ) 2 when the y-increment 
is changed. An expression for determining r such that M remains within 
reasonable bounds can be derived by relating the value of N, N , at which 
this change occurs to the value of N which gives the desired value of 
N » Equation (33) is 


An. 


M x (r t - 1) 


N r 


- 1 


V4 


Also assume 
M 


= C 


N 


(42) 


1, N 


where C is a constant. The most logical choice for this constant is 
probably 4. Since r^ c - r^ c 1 , Eqs. (33) and (42) give: 


r 

t 


(An r(c - i) 

c 


+ 1 


(43) 


If this value of r^ is used* M 1 N will vary between and M^/C for 
m = 1. The asymptotic , minimum M 1 N can easily be shown to be 


lira 
m •* 00 


1, N 



(44) 


Due to the variation of r with m, I is growing progressively smaller* 

\ y 

A slight modification which definitely warrants study would be to keep 
I fixed. Every other grid point would then be deleted based on r\ ; in- 
stead of Ar) . Initially, the two procedures would give identical re- 
sults; however, as m increases significant differences in the various 
parameters would arise. Smaller values of Ap™ would result in order to 

keep I constant as r becomes progressively larger. 

y 
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The progression of r y and I for two different values of Art c are 
given in Table 2.2. The error curves which were obtained with the parame- 
ters given in Table 2.2 for the problem of a semi-infinite solid with 
a step change in surface temperature are given in Figs. 2.5 and 2.6. 

Since the error curves for these cases lie so close together even with 
the greatly expanded scale used, the spread of the seven additional 
curves between An and approximately 2Ap are simply indicated by a 

C C 

cross-hatched area. The error was always less than 0.9 percent and 
greater than 0.2 percent with Ari c equal to 0.15 and for N between 16 
and 289, The maximum and minimum errors were reduced to 0,4 percent 
and 0,1 percent with Ari c = 0.10 and N between 37 and 369, M varied be- 
tween 2,0 and 0.5 for the first error curve (m = 1), and it varied be- 
tween 2,8 and 0,67 (see Eq, (44) ) for the last error curve 9 m = 7, 

The corresponding range of for the two cases is approximately from 
21 to 4 X 10 5 and 48 to 9 X 10 5 . A uniform grid with the initial 
spatial and time increments would have required approximately 200 3 000 
times more computational effort! 

Figure 2,7 shows a set of error curves for the Blasius problem. 

The parameters used are the same as those employed in Case I of the 
conduction problem and Table 2,2 is applicable. It is noted that the 
larger values of r^ have a much smaller influence on the accuracy if 
the calculation is begun with small values of r^. For example 9 compare 
the results of Figs, 2,1 and 2,2 with Figs, 2,5 through 2.7, 

Hie tremendous savings in computational effort makes the general 
procedure look indeed attractive. Much further study is required and 
general conclusions cannot be drawn at this time. It is clear that the 
possible savings in computational effort in the solution of turbulent 
boundary layer problems and other nonsimilar problems is not merely a 
factor of 2 or 3 but several orders of magnitude. 


28 



















Si 



0 <><>\ 




'<X, 

V 


\\ 


J 


c** 

<1 


C\'1 

d 


fl ^ A ^ : 




| iZ'i Jiij 

“ 51 !1 o 


fl ] p ^ 


0 


Figure 2.5 The Influence of Increment Changes on Accuracy 3 Ari c = 0.15 (Semi- 
Infinite Solid with Step Change in Surface Temperature) 
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Figure 2.7 The Influence of Increment Change on Accuracy 9 An =0.15 (The Blasius Problem) 



3. FINITE DIFFERENCE ALGORITHMS WITH REPRESENTATIVE RESULTS 


3.1 Algorithms for Compressible Flows of the Boundary Layer Type 

As a result of the comprehensive study of the application of finite 
difference methods , three different algorithms were developed and employed 
in numerical experiments. These algorithms were presented in detail in an 
earlier report [1]. Many solutions of laminar incompressible flows were 
also given in [1], and pertinent results in graphical and tabular forms 
were presented in order to establish the accuracy of the algorithms. 

The advantages and disadvantages of each algorithm were discussed and, 
based on these considerations, the Crank-Nicolson algorithm was selected 
to solve the film cooling problem. Although the accuracy of the Crank- 
Nicolson method for compressible flows was demonstrated in a status re- 
port, the difference equations and the associated solution procedures were 
not presented. Hence, they will be given in this section. 


The Crank-Nicolson method was employed to solve the film cooling 
problem for two reasons: First, this implicit method eliminates the prob- 

lem of instability associated with explicit methods. Although the scheme 
results in simultaneous nonlinear equations, their solution does not generally 
require more computational effort per step than the explicit method. Second, 
the average differences generally result in greater accuracy compared to the 
Standard Implicit and Explicit methods. Some representative results showing 
the accuracy of the Crank-Nicolson algorithm for the solution of laminar 
incompressible and compressible boundary layer flows are given in Figs. 4 
and 5 of [3], 


The solution procedure associated with the Crank-Nicolson method w ill 
be described by an example. Consider the case of compressible boundary layer 
flow over a flat plate. Dimensionless variables are introduced in order to 
leave u^, T^, p^, and k^ unspecified. Specifically, 




oo 


T* 



, and k* 

•no 



The governing difference equations derived from a physical model (see [1]) 
in terms of the dimensionless variables are: 
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Continuity: 


( P* V ' f \,n*» = (pVt) i-l,n +14 


( P* U *)i-l,n + ( P V) i,n - ( P* u&) i _ x> 


n+1 


( p*u* ) 


n + 1 


(45) 


Momentum: 


(o*u*) (u* - un ) 

i, n+Vi i , n+1 i f n 


+ (p*v» t ). 


i n+Vz 


U -U „ + U. , -U, 

i + l f n+1 i -1, n+1 i + 1 , n l -1 1 n 


1 

2M 


[*. 


'A, n+ 1 ^ U iV 1 , n+ 1 U f, n+ 1 ^ + W'-A, n+ 1 ^ U i ' - 1 , 


- u* ) 

n+1 1 , n+ 1 


+ y* „ (u* , - u.“ ) + y* (u*- , - u*. ) 

i + Vz, n i + l, n i , n 1 - Vz , n 1 - 1 , n i , n 


(46) 


Energy : 


(p*u*). XI/ (T* - T* ) + (p-V 1 ). ' 

^ i , n + Vi i , n+ 1 i , n i , n+Vz 


T* - T* + T* - T* 

i + l, n+1 i-l,n+l i+l,n 1-1, 


I 

I 

* 


2(Pr)(M) 


V* ( T _ T * * ) 

i+H n+1 i + 1, n+1 i , n + 1 ' 


. i a ( t; t : _ T;’; ) + V* (T* - T*' ; ) 

+ -V h n+1 ^ l -1, n+ 1 i, n.+ 1 ' K i + l A t n i + l, n i , n 


+ k* „ (T* - T. 

i -‘A, n i -1, n 


T, n >] 


( T - . J-L Mj-V (U* 

16 (M) n , n + 'A i + l, n+1 


Equation of State : 
p*T* = 1 

Property Relationship : 


u* xl + u* - u* ) (47 )_ 

1 - 1 , n+1 i + l, n l -1, n H 


(48) 


' y* = k* = (T*) w 


(49) 


where Pr and Ma are the Prandtl and Mach numbers, respectively. The quanti- 
ties (v 1 )^ = (v/u ) (Ax/Ay), and M = [(Ay) 2 u co ]/(v ro Ax) were introduced in 
order that the spatial increments Ax and Ay could remain unspecified. The 
indices i ,n represent the location y = (i - l)Ay and x n = n(Ax). Since u 
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and T are defined only at the discrete locations y, and x q 3 the quantities 

which must be known at the locations y. .. or x ,, are calculated as follows: 

-'i +Vi n + Vz 

1 


(P¥| I,..» 3 2 C(,¥1 I,.„ + 
u* = k* = ( . 

+ Vi, n i + Vzj n l 


T* , + T* 

i + 1 , n l f n 


w* 


T* , + T* 

i , n+ X l , n 


"i, n + Vz \ 2 

The initial and boundary conditions are 


u* =1 and T* = 1 , i > 0 

i,o i ; o 


u* = v ! =0 and T* = 1 , n > 0 

o, n o, n o.n 5 — 


u* = 1 and T* = 1 , n > 0 , i = I (where I is large) 

1 » n i , n 

An examination of Eqs. (46) and (47) shows that if the quantities 

( P* U * ) i,n + %» ( P* V ' ,) i I n + V4» y f + -A, n + 1 > md k *+>A, n+ 1 Were knOWn » the equations 
would be linear. If these "unknown" quantities were approximated, Eqs. (46) 

and (47) could be solved iteratively, and the results from the iteration 
could be employed to correct these "unknowns " after each iteration. Numeri- 
cal experiments indicated that the iteration process converges rapidly if 
the initial approximation to these "unknowns" is accurate. 

A predictor for u* , and T* , was derived in [1] and is given by 

l , n+l i, n+1 ° J 

<t>° .. = 3$. - 3(f). + <j). , n > 1 (50a) 

where <J) denotes either u* or T*. The expression is valid only for n > 1 
and approximations for u* and T* are also required for n equal to zero and 
one. A linear extrapolation is used for n = 1; that is 


6 = 26 - 6 

i , 2 y i , 1 U , o 

and the initial condition is used for n = 0; that is 


(50b) 


<f) = 4> 

i i o 


(50c ) 
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The approximation for <J>. is arbitrary so long as the iteration converges. 
Equation (50a) is identical to that which is obtained by equating the 
second central differences 

6 2 <j). = 6 2 4>. 

x i , n x i , n -1 

2 2 

Thus , if 6 d>. and 5 <f>. , are of opposite sign, Eq. (50a) would give 

3 x i, n x i , n -1 ° 

an inaccurate prediction for (J). n+1 * This difficulty had been encountered 
in some laminar and turbulent problems . Several other predictors were 
tested and two were found which avoided this difficulty and also gave ac- 
curate predictions. One of them is the Standard Explicit method (see [1]) 
but its application is limited to incompressible problems. The other one is 


6 2 <J>. = S 2 <f>. 9 

x i , n x T i , n -2 


or 


d>° = d). +6 (j) w + 6 d>. - 6 c n , n > 2 (51) 

Y i , n+ 1 Y i , n x Y i , n -Vz x Y i , n -3/ 2 x Y i , n -5 / 2 3 


For n <_ 2 , 6^ <}>. , <$„ <fc 


and 6 d>. „ , „ were assumed to be zero, 

i , -3/ 2 x Y i , -5/ 2 


These assumptions are arbitrary as long as the solution converges. 
Equation (51) does not require additional computational effort compared 
to Eq. (50a) if the central differences are stored. 


A method of solution which is similar to the predictor-corrector me- 
thods employed in the solution of ordinary differential equations is used. 

A "predictor" , Eqs. (50) or (51), is employed to predict u* n+1 and T* n+ 1 * 
Next 9 the continuity equation , the equation of state, and the property re- 
lationships are used to obtain (p*v** ). n+1/ ? p* , y* , and k*. Then the- mo- 
mentum and energy equations, the "correctors", are used iteratively to 
correct the discrete velocities and temperatures, and to add the required 

stability to the overall process. The continuity equation is employed 

t 

after each application of either "corrector" to correct (p*v* ). n+1/ . y, 

p, and k are recalculated only after each application of the energy 
equation. Since the predictor provides accurate estimates of the required 
dependent variables after the first few steps, the correctors usually must 
be applied only once. 

After the quantities p*u*, p*v** , y*, and k* are calculated, Eqs. (46) 
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and (47) can be rewritten in the form 


A ± K n .K „ X K 

A. 6 + B d) +. C cb 

i M - 1 , n + 1 i M , n+ 1 i M + l, 


A = b. , 1 < i < I 

n+ 1 i 5 . — 


( 52 ) 


where A , B i , C. , and b are all known quantities. The superscript K de- 
notes the kth approximation obtained from the corrector; K = 0 corresponds 
to the prediction. Thus, the momentum and the energy equations are re- 
duced to two sets of linear algebraic equations with tridiagonal co- 
efficient matrices. The simultaneous equations are solved by employing 
the algorithm of Section 3.1 of [6]. 


Since the corrector must be applied iteratively, a means of terminat- 
ing this iterative process is required. A -test obtained from the con- 
servation of mass is 


(100 percent) 


K 

r (p-iv* > f 

L M ; I , n + Vzj 


K-l 


< e 


(53a) 


The expression in the left-hand side of Eq. (53a) is equal to 1/Ay times 
the change between iterations of mass flow out of the control volume across 
the surface perpendicular to the plate at location n + 1. To reduce the 
probability of accidentally satisfying Criterion (53a), a second criterion 


(100 percent) 


. (p * U * ) 2,n + l] K ' [ ( P i ’ !U * ) 2 I n + l] K ’ 1 


[ ( P* U;,:) 2,n + l_ 


K 


< o.i e 


(53b) 


is introduced. Criterion (53b) is generally less stringent than Criterion 
(53a). Numerical results showed that £ = 1 is sufficiently stringent for 
engineering calculations. 


Finally, some of the advantages of the algorithms being developed in 
this study in contrast to many of those described in the literature are: 
1. The solution is effected in the physical plane. This makes 
the algorithms more general and enables their adaptation to a 
variety of problems with little or no modification. The physi- 
cal variables are available for direct examination; hence, the 
characteristics of the problem are more easily discovered. 
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No initial condition is required in v or 9u/9x. Many tech- 
niques reported in the literature require such conditions. 
These conditions are incompatible with the boundary layer 
approximations. 

3. The finite difference equations are based on a physical 
model; hence, the derivation of expressions for the wall 
shear stress, the wall heat flux, and the boundary conditions 
can be accomplished in a rational manner which is consistent 
with the approximations employed in the interior of the flow 
field. 


3.2 A Turbulent Model and the Associated Finite Difference Procedure 

Consider the case of two-dimensional, incompressible, turbulent flow, 
the governing mean flow equations employing an eddy viscosity concept 
are (see [4]): 


9u 9v 
9x 9y 


0 


(54) 


(— 9u — 9u \ dp 9f/ A v9ul 

P ( U 97 3y j“ dx 3y [ (p V 8y J (55) 

where u, v are the mean velocities, u = u + u' , p is the mean pressure, and 

is the turbulent mixing coefficient which is often called the "eddy 11 vis 

cosity. A is defined as 
J T 


A = ~P U ' yT 
T du/dy 

where u 1 , v ! are the velocities of fluctuations. Equations (54) and (55) 
are made to correspond to the incompressible, laminar boundary layer 
equations by introducing an effective viscosity 

h e = U + \ (56) 


Hence, the finite difference algorithms of [l] can be employed to effect 
a solution. A is not a property of the fluid like y, but is dependent 
upon the mean velocity u. In this preliminary investigation, Prandtl f s 
mixing length hypothesis was employed in the viscosity model, thus, Eq. 
(56 ) becomes 
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(57) 


V 


e 


= y + p& 2 


3u 

3y 


where i is the mixing length. The mixing length expressions employed are 
given by: 


= 0.41 [1 - exp (~y + /26)] (y/<$), y/6 £ 0.1 (58a) 

[0.41 (y/6) - 1.53506 C 2 + 2.75625 ; 3 - 1.88425 £ 4 ] 

[1 - exp (-y + /26 ) , 0,1 £ y/6 £0.6 (58b) 


l 


6 


0.089 , 


y/6 >0.6 


(58c) 


where 



C = (y/6 - 0.1) 


and 6 is the boundary layer thickness. Equation (58) is a modification of 
Piet cherts expressions [7]. Pletcher T s expressions were modified in order 
to avoid the discontinuity in Z/6 at y/6 = 0.1 which occurs if y + is small. 


If the turbulent boundary layer equations were to be solved by a 
finite difference technique employing a uniform network, small increments 
are required in order to include grid points in the laminar sub-layer. How- 
ever, the thickness of this layer is small relative to the thickness of 

7 + 

the boundary layer; for Re^ = 1.87 X 10 , = 12 ,000 [4], if an in- 

crement of Ay + = 4 is employed, it would require 3,000 grid points in the 
normal direction. Hence, a nonuniform network must be used. This con- 
sideration provided the impetus for the study presented in Section 2.3. 

The finite difference equations and the associated solution procedure 
for the case of turbulent flow over a flat plate will be considered. 

The nohuniform increments are given by the following geometrical progres- 
sions. 

Ay. = r 1 Ay, , i > 1 (22) 

i y 1 — 

and 

Ax = r” 1 Ax , n > 1 (59) 

n x 1 — 
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Therefore , 


^ i [ "V~~ ) • 


i > 1 


(24) 


and 


/r" - 1 

V AX - 


n > 0 


(60) 


The governing difference equations obtained by employing the Crank - 
Nicolson method are: 

Continuity: 


< V ' >i\ n+V4 


_ >i - 1 , 


n+Vi ^ 1 

r + 2( u f-i,« 

y 


+ u* - u* - u?’ c ) (61) 

i f n i -1 f n + X i f n+ 1 / 


Momentum: 


u? 


(v' ) f 

(S* - u* ) + 


i , n + V4 v “i,n+l “i , n 7 ' 2r^ (1 + r y ) ^L^i+l.n+l y 


2 2 

+ (r - 1) u* , - r u* 

i, n+l y i 


" 2 2 

u* + (r - 1) u* - r u; 
» + 1 . n y • - " • 


i , n y 




(1 + r )(r 2 ) 2 M 

y y i 


W 


i + ‘/2, n+ 1 


( U* , , - U?’ ! , ) + U* ( U* , 

' l+l, n+l » i n+ 1 ' l - l A, + 1 v i-l, 


- U* ) 
n+l l , n+ 1 1 


u« 

( " U.° ) + yb ( U« , - Ub ) 

r ' 1 + 1 , n i , n • i n v 1 -1 , n 1 , n ' 


(62) 


where 


u ;; = — , (v* )J 


u 

u 


V i,n + >/ 2 / AX n+l 


i , n + !/2 


U 


Ay i -i, 


n+ 1 


pu 

, ^oo 0 

y* = l + l 2 

y 


3u*“ 


3y 


and 


u J A y! ) 

M E 


1 v(Ax x ) 


( 63 ) 


-1, n+ 1 
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Lastly, the expressions for y* and £/6 in difference form are: 

^.*1 = 1 + (Re '. v* (?) 2 

* i+^n+l 


(i 6 - i) 


u ?\ . - ir 

i + 1 f n+ 1 i , n 


i -X 


(64) 


(i). 


i + n + 1 


0.41 [1 - exp (-y. +1/2i n+l /26)] 


(!). 


i +Vi, n + 1 

< 0.1 ' 


(65a) 


i + V 2 , n+ 1 


(I). „ ,= I" 0 - 41 !?). „ , - 1 - 53506 ?U».i 2 ' 75625 


i +Vz t n+ 1 


i +Vzj n + X 


i + Vzj n+ 1 


1.88425 C 


i +V 2 , n+ X 


1 - exp <-y^ n+1 /26) 


where 


(?) = °-° 89 . 

i +Vz, n+1 


0.1 < (y/S ). +l/ji n+1 10.6 
> 0.6 


y 

J 


i +Vz, n+ 1 


t U 00 

Re i = V 


(65b) 

(65c) 


(66) 


i -X 


W + i = ( Re4 *>* £ i- 


/y\ 

UA 


X X y 1 , n+ X 


[i - 1 + (r _1 /2)] 
e y 


- 1 + 


i +V 2 , n+ 1 


jT r~ iy 


4*K,.+ 1 =[(?) “ 0 - 1 ] 

t *- 1 + t 2| n+ X -J 


r 1 - l v 
y 


i = 1 +| , 

e \ r - 1 


(26) 


ig is the location corresponding to the edge of the boundary layer, it is 
usually found by interpolation because the velocities are known only at 
discrete locations. is the operator employed in calculating the 


41 



derivatives at the wall; it is given by* 

C 4 *4 + C 3 ^3 + C 2 ^2 + C 1 

A y *, (f 

where 

C = -r 3 (r + 2)(2r 2 + 4r +3) 

1 y y y y 

6 5 4 3 2 

C = 2r + 8r + 13r + Hr + 5r + 2r 

2 y y y y y y 

C = -(2r* + 5r 3 + 5r 2 + 3r +2) 

3 y y y y 

C = r +2 

4 y 

C = 2r 3 (r + l) 3 

5 y y 


( 67 ) 


The initial and boundary conditions are the same as that of the Blasius 
problem and are not repeated. 

An examination of the difference equations show that neither Ax nor 
Ay needs to be specified to effect the solution; 9 Re^ 5 r^ 5 r y 9 I, and 
N (or their equivalent values 9 Eqs . (26) and (27) ) are the only parame- 
ters needed. The solution procedure employed is the same as that de- 
scribed in [1] with the following exceptions: 

1. The "predictor" used is Eq. (51) instead of Eq. (50). Although 
Eq. (51) was derived for uniform x- increment,, it provided ac- 
curate predictions for u. n+1 for the nonuniform x-increment 
cases as expected. 

2. The effective viscosity must be recalculated after the appli- 
cation of the "corrector" 9 because the eddy viscosity is a 
function of u. 


Figure 3.1 shows the influence of Re^ 9 r^ 5 and r^ on the accuracy 
of the local skin friction coefficient 9 c f ! . The results are preliminary 
further investigation of the influence of M and increment changes is 
planned. 


*See [1] for the derivation of this operator for r^ = 1. 
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Figure 3.1 The Influence of r^ , r y and on the Accuracy of 
Incompressible Turbulent Boundary Layer Flow 
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4. RESULTS FROM THE STUDY OF FILM COOLING 

Results for laminar incompressible film cooling neglecting viscous 
dissipation were presented in [1]. Results obtained with a simplified 
model for laminar gaseous film cooling including the influences of 
compressibility and viscous dissipation will be presented here. These 
results are presented in generalized form, where the generalized in- 
dependent variables were obtained by the techniques given in [3]. 

Gaseous film cooling involves the mixing of two moving streams 
in the vicinity of a solid wall. A coolant or secondary stream with 
a velocity u g and a temperature is ejected through a slot in a 
direction tangent to the surface. The objective is to protect the 
surface from the high temperature external stream, the primary stream, 
whose velocity and temperature are u^ and T^ . The primary flow in 
the case of interest is a supersonic and semi-infinite stream. It 
is assumed that a converging nozzle is employed for the coolant flow; 
hence, the velocity at the exit of the slot is either sonic or sub- 
sonic. 

Figure 4.1 shows the flow configuration studied. This flow con- 
figuration results when the pressure of the secondary stream at the 
slot exit is equal to the pressure of the primary stream. The split- 
ter plate is assumed to have zero thickness; hence, no expansion fans 
or shock waves will be present. The following assumptions were em- 
ployed; 

* 1. The flow field is laminar, two-dimensional and steady. 

2. The secondary and primary fluids are the same and are 
ideal fluids. All relevant fluid properties are functions 
of temperature only. 

3. At the exit of the slot, the secondary flow is either sub- 
sonic or sonic and the primary flow is supersonic. The 
static pressure of the secondary flow at this location is 
equal to the static pressure of the primary flow; hence, 
the flow configuration of Fig. 4.1 is obtained. 

4. The boundary layer equations are valid in both regions I 
and II, (These regions are defined in Fig. 4.1.) 
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5. The free stream pressure and temperature do not vary 
in a direction along the film cooled surface. 

6. The initial velocity and temperature profiles of both 
streams are uniform and the velocity vectors are 
parallel. 

7. The splitter plate which initially separates the two 
streams has zero thickness at x = 0. 

8. Heat transfer by thermal radiation is negligible. 

9. The surface is isothermal and at the temperature of 
the secondary stream, T g . 

The governing equations with these assumptions are: 

Continuity: 


(pu) + |^r (pv) = 0 


Momentum: 



Equation of State: 


pT = 2- = a constant 
K 

Property Relations : 

y = y(T) , k = k(T) , and c p = c p (T) 
The initial boundary conditions are 


at x = 0 , 

u = u and 

s 

T = 

T s> 

y 

< s 


u = u and 

P 

T = 

V 

y 

> s 


at y = 0 u=v=0 and T = T 

S 

x > 0 

at y + “ u=u = u and T = T = T 

OO P OO p 


( 68 ) 
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The property relations remain to be specified. It will be as- 
sumed that the Prandti number and the specific heat are constants. 

Thus , y is proportional to k and a single function needs to be speci- 
fied. The viscosity can be accurately specified by Sutherland's formula 



(H 


3/ 2 T + S 
o I 


T + S 


(73) 


where ]i is the viscosity at the reference temperature T which can be 

o o 

chosen arbitrarily and is Sutherland's constant. On the other hand, 
a simple power law 



(74) 


is often used because this less accurate law gives rise to one less 
parameter if T/T q is employed as the dimensionless temperature. For 
this reason it was used in this investigation. A more accurate ex- 
pression $ such as Sutherland 1 s formula , as well as variable specific 
heat could be easily incorporated into the finite difference algorithm. 


Figures 4.2 through 4.5 show the results of compressible film 
cooling. In order to illustrate vividly the effectiveness of cooling , 
the shear stress , the heat flux at the wall, and the total rate of 
heat transfer are normalized by the corresponding values for zero slot 
height . Specifically , 


T 

_ w 


a - % 


~ 

w 

n ’ ^ <c 

5 q " q 



s=0 

S=0 

S=0 


The results in these figures were obtained with a) = 1, y = 1.4, and 
Pr = 1. These values were chosen only because x | , q"| s _Q, an( ^ 

q| s= () are known ^ or "this case and do not have to be obtained numeri- 
cally. 


Figure 4.2 shows a family of temperature profiles. The parameter 
is the dimensionless x-location, £ = (x/s Hv^/u s ) . The change in the 
heating rate can be clearly seen by examing the slopes at y/s equal to 
zero. Figure 4.3 shows the variation of the normalized heat flux at 
the wall with the dimensionless x-location , £. The velocity ratio. 
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/u p , is 0.25 and the Mach number of the primary stream, Ma p , is 
equal to two; the parameter is the temperature ratio, T s /T . 

Figure 4.4 shows the same results including, in this case, also the 
the normalized shear stress and heat flow rate with a new dimension- 
less x-coordinate , defined as £ = (x/S)(v /u S) = (T /T ) 

The normalized shear stress, heat flux, and total heat flow rate be- 
come, for practical purposes, independent of T s / T p when E >1 is used 
for the abscissa in place of 

Figure 4.5 shows the variation of the normalized wall shear 
stress, wall heat flux, and total heat flow rate with the dimension- 
less x-location, = (x/S) (v /uJB). The velocity ratio is 0.2 and 
the temperature ratio is 0.3; the parameter is the Mach number. The 
plots show the influence of Mach number on the heat flow ratio to 
be small. 

In this analysis, the primary stream is assumed to be uniform 
which is unrealistic in most applications. Hence, the results are 
upper limits of the more practical case with an initial boundary 
layer in the primary stream. Furthermore, the assumptions of laminar 
flow is quite restrictive. The mixing region, in most practical cases 
will be turbulent. The boundary layer which forms along the wall will 
be laminar inti ally but will become turbulent at a relatively low 
Reynolds number due to the turbulence in the neighboring mixing region 
Hence, the turbulent case is now being investigated. The preliminary 
results of Sections 2.3 and 3.2 stemmed from this effort. Some of the 
other assumptions currently being used also need to be examined more 
closely. 
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